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The energy reconstruction of extensive air showers measured with the LOEAR Radboud Air Shower Array (LORA) is presented in 
detail. LORA is a particle detector array located in the center of the LOEAR radio telescope in the Netherlands. The aim of this 
work is to provide an accurate and independent energy measurement for the air showers measured through their radio signal with 
the LOEAR antennas. The energy reconstruction is performed using a parameterized relation between the measured shower size 
and the cosmic-ray energy obtained from air shower simulations. In order to illustrate the capabilities of LORA, the all-particle 
cosmic-ray energy spectrum has been reconstructed, assuming that cosmic rays are composed only of protons or iron nuclei in 
the energy range between ^ 2 x 10^® and 2 x 10^® eV. The results are compatible with literature values and a changing mass 
composition in the transition region from a galactic to an extragalactic origin of cosmic rays. 
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1. Introduction 


The quest for the origin of cosmic rays is one of the most fun¬ 
damental problems in Astroparticle Physics lillS] . Since the 
discovery of these highly energetic particles more than a cen¬ 
tury ago, numerous measurements of several of their properties 
have been made, using sophisticated instruments (see e.g. Ref. 
H for a review). However, the exact nature of their sources still 
remains an open question. The search is mainly hindered due to 
the fact that cosmic rays, being charged particles, are scattered 
or deflected by the Galactic and inter-galactic magnetic fields 
during their propagation to the Earth, making it extremely diffi¬ 
cult to reconstruct the direction of their sources. Nevertheless, 
observed cosmic-ray properties like the energy spectrum and 
composition have been used to understand and characterize the 
properties of the sources such as their Galactic or extragalactic 
nature, the cosmic-ray production spectrum and the power in¬ 
jected into cosmic rays (see e.g. Refs. 10,0, S SS[13 [HI for 
recent reviews). 

LOEAR, the LOw Erequency ARray, is an astronomical radio 
telescope 113] . It has been designed to measure the properties of 
cosmic rays above ^ 10^® eV by detecting radio emission from 
extensive air showers in the frequency range of 10 — 240 MHz 
1 13r . One of the main goals of the LOEAR key science project 
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Cosmic Rays is to provide an accurate measurement of the mass 
composition of cosmic rays in the energy range between ^10^® 
and ^ 10^® eV, a region where the transition from Galactic to 
extragalactic cosmic rays is expected. This is being carried out 
by measuring the depth of the shower maximum (Ai„iax), using 
a technique based on the reconstruction of the two-dimensional 
radio intensity profile on the ground QQ. Another focus 
of the LOEAR cosmic-ray measurements is to understand the 
nature and production mechanisms of the radio emission from 
air showers. This is done by measuring various properties of the 
radio signals in great detail such as their polarization properties, 
the radio wave front and relativistic time compression effects on 
the emission profile 


In order to assist the radio measurement of air showers with 
LOEAR, we have built a particle detector array LORA (LOEAR 
Radboud Air Shower Array) in the center of LOEAR in . Its 
main objectives are to trigger the read-out of the LOEAR ra¬ 
dio antennas to register radio signals from air showers, and to 
provide basic air shower parameters such as the position of the 
shower axis as well as the energy and the arrival direction of 
the incoming cosmic-ray. These parameters are used to cross¬ 
check the reconstruction of air shower properties, based on the 
measured radio signals. Currently, given the lack of an absolute 
calibration of the radio signals, the cosmic-ray energy is esti¬ 
mated through the reconstruction of the particle data. There¬ 
fore, an accurate energy reconstruction with LORA is essential 
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Figure 1: Layout of the LORA array in the LOFAR core. The filled black 
squares represent the LORA detectors, the crosses the LOFAR low-band an¬ 
tennas and the empty squares the high-band antennas. The dashed circle in the 
figure illustrates the fiducial area of a radius of 150 m, which is used in the 
analysis. 


for a proper understanding of the air showers measured with 
LOFAR. 

In this article, we describe in detail the various steps of the 
energy reconstruction and present the cosmic-ray energy spec¬ 
trum above ^ 10^® eV as measured with LORA. The article is 
organized as follows. A short description of the set-up will be 
given in Section |2] followed by a description of the data anal¬ 
ysis technique in Section [3] The various steps involved in the 
Monte-Carlo simulation studies of the array will be described 
in SectionlH and a comparison between measurements and sim¬ 
ulations for some of the air shower properties will be given 
in Section |5] In Sections and [T] the energy calibration, the 
uncertainties in the reconstructed energies, and reconstructed 
cosmic-ray intensity will be described. The measured cosmic- 
ray spectrum and a comparison with the measurements of other 
experiments will be presented in Section|8] followed by a short 
conclusion and a future outlook. 


2. LORA experimental set-up and operation 


LORA (the LOFAR Radboud Air Shower Array) consists 
of an array of 20 plastic scintillation detectors of size 
0.95 m X 0.95 m each, distributed over a circular area with 
a diameter of ~ 320 m in the center of LOFAR lIlQIl . The array 
is subdivided into 5 units, each comprising of 4 detectors. The 
detectors have a spacing between 50 — 100 m, and have been de¬ 
signed to measure cosmic rays with energies above ^ 10^® eV. 


The array is co-located with six LOFAR station^ The layout 
of the array is shown in Figure [1] The data acquisition in each 
unit is controlled locally. A local trigger condition of 3 out of 
4 detectors is set for each unit, and an event is accepted for a 
read-out of the full array when at least one unit has been trig¬ 
gered. A high-level trigger for the LOFAR radio antennas is 
formed when at least 13 out of the 20 detectors have measured 
a signal above threshold. More technical details can be found 
in Ref. d. 


3. Data selection and analysis 

Data collected with the LORA array since its first science 
operation in June 2011 until October 2014 are used. Only data 
collected in periods with all 20 detectors in operation will be 
considered. This amounts to a total of 706.9 days of data. For 
the analysis, only showers that trigger a minimum of 5 detectors 
will be considered, which corresponds to a total of 1,861, 045 
air showers. 

For every measured shower, the signal arrival time and the 
energy deposit in each detector are recorded. The relative sig¬ 
nal arrival times between the detectors are used to reconstruct 
the arrival direction of the primary cosmic ray. The energy de¬ 
posits are used to reconstruct the position of the shower axis 
and the shower size (the effective number of charged particles 
at the ground). The latter is determined in terms of the number 
of vertical equivalent charged particles, which may also include 
converted photons in addition to the dominant charged particles 
- electrons and muons. The shower axis position and the shower 
size are determined simultaneously by fitting a lateral density 
distribution function to the measured two-dimensional distribu¬ 
tion of particle densities, projected into the shower plane. The 
particle density in each detector is obtained by first dividing the 
track-length-correctec0 energy deposition by the energy depo¬ 
sition of a single particle obtained from calibration, and then 
by further dividing by the projected area of the detector in the 
shower plane. The lateral density distribution of an air shower is 
generally described by the Nishimura-Kamata-Greisen (NKG) 
function which is given by 21 ] 


p(r) = NchC{s) ( — 
\rM 


s-2 


1 + 


tm) 


s—4.5 


( 1 ) 


where p(r) represents the particle density in the shower plane at 
a radial distance r from the shower axis, Nch is the shower size, 
s is shower age or lateral shape parameter and tm is the radius 
parameter which is basically a measure of the lateral spread of 
the shower. The function C(s) is given by 


C{s) = 


r(4.5-s) 


27rr^r(s)r(4.5-2s)- 


( 2 ) 


'Each LOFAR station consists of 96 low-band and 48 high-band antennas, 
operating in the frequency range of 10 — 80 MHz and 110 — 240 MHz respec¬ 
tively. 

^The measured energy deposit in each detector is corrected for the increase 
in the path length of the incident particles through the detector by multiply¬ 
ing by a cos 6 factor where 9 is the zenith angle of the reconstructed arrival 
direction of the primary cosmic ray. 
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Figure 2: Left: Normalized distribution of radius pai'ameter, for the measured showers with reconstmcted size logj^Q Nch > 6.40 and zenith angles in the 
range 0° — 15°. The inset shows a Gaussian fit (represented by the line) to the distribution around the maximum. Right: Averaged lateral distribution of measured 
showers with a reconstructed size in the range of 6.40 < logj^Q A^ch < 7.75. Only uncertainties for the uppermost and the lower-most distribution are shown. The 
lines represent fits of an NKG function, keeping the shower age parameter fixed at 5 = 1.7. 


In the case of LORA, the value of tm is determined from the fit 
along with A^ch and the position of the shower axis. The param¬ 
eter s is kept constant at a value of 1.7 throughout the fitting 
process. Simultaneous fitting of both tm and s results in fits 
of poorer quality. Simulation studies have shown that keeping s 
constant gives better results than keeping tm constant 1221] . The 
fitting procedure is repeated three times with the output of each 
fit taken as starting values for the next iteration. Details about 
the minimization procedure and the choice of starting values, as 
well as the reconstruction of the arrival direction of the primary 
particle are described in Ref. S. 

All showers that trigger at least 5 detectors with a minimum 
of 1 particle m~^ are allowed to pass through the reconstruction 
algorithm, and their shower parameters are calculated. Further¬ 
more, only showers whose reconstructed position of the shower 
axis falls within 150 m from the center of the array are se¬ 
lected. The normalized distribution of tm values for the se¬ 
lected showers with reconstructed sizes logj^g -^ch > 6.40 and 
reconstructed zenith angles in the range of 0° — 15° are shown 
in Figure |2] (left panel). The inset shows a closer view of the 
distribution around the maximum value between 12 and 48 m, 
and a Gaussian fit to the distribution. The fit gives a peak value 
of tm = 30.33 ± 0.13 m. Figure |2] (right panel) shows the aver¬ 
aged lateral distributions of the measured showers for different 
reconstructed size bins in the range of 6.40 < logj^Q Vch < 7-75 
for zenith angle between 0° and 15°. The distributions include 
events that passed through the same selection cuts applied in 
the left panel of Figure |2] and have tm values in the range of 
10—200 m. The averaged distributions are obtained by stacking 
together the lateral distributions of all individual showers con¬ 
tained in each size bin. The lines in the plot represent the fits to 
the data using ([T]). To avoid clumsiness of the plots, uncertain¬ 
ties are shown only for the size bins of logj^p ^ch = 6.40 — 6.55 
and log^o -Vch = 7.60 — 7.75. However, all the respective un¬ 
certainties are taken into account in the fitting procedure. The 
values of tm obtained from the fits are in the range of 23 — 31 m. 

The shower size gives a good measure of the energy of the 


primary cosmic-ray particle, initiating the air shower. There¬ 
fore, the shower size distribution should reflect the energy dis¬ 
tribution of the cosmic rays at size values where the primary en¬ 
ergy is above the detector threshold. Figure [3] shows the distri¬ 
bution of reconstructed shower sizes for all the measured show¬ 
ers that passed through the various trigger and quality cuts ap¬ 
plied in the analysis. This corresponds to a total of 322, 664 air 
showers. The distribution shows a steep rise as Nch increases 
which is due to the sharp increase in the detector acceptance 
(see section m as function of the primary energy. After reach¬ 
ing a maximum, the distribution falls off steeply which is due 
to the power-law behavior of the cosmic-ray spectrum. The 
peak of the total distribution gives the shower size threshold of 
the detector array. Fitting a Gaussian function around the peak 
gives a value of logj^g A^ch = 5.92. Also shown in Figure[2are 
the reconstructed size distributions for four zenith angle bins: 
0° - 15°, 15° - 24°, 24° - 30°, and 30° - 35°. The parameter¬ 
ization of the cosmic-ray energy will be determined separately 
for each zenith angle bin (see Section |6|. All cuts applied in 
this analysis are summarized in Table[T] 


4. Simulations 

Detailed simulation studies have been carried out in order to 
understand the performance of the array and to determine var¬ 
ious characteristics of the array, such as the trigger and recon¬ 
struction efficiencies, the reconstruction accuracies of shower 
parameters, the relation between reconstructed size and primary 
energy, and the accuracy in the energy reconstruction. In this 
section, the various steps involved in the simulations will be 
described. 


4.1. Air shower simulations 


Air showers are simulated using the CORSIKA simulation 
package (version 7.4387) 1231] . The interactions of hadronic 
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Table 1: Selection cuts applied in the analysis of both the measurements and the air shower simulations. 


Trigger condition: 

Single unit trigger: 

Analysis: 

Number of leftover showers: 

3/4 detectors 

5 detectors with > 1 particle m~^ 
1,861,045 

Quality cuts: 


Zenith angle: 

e < 35° 

Position of the shower axis: 

< 150 m from array center 

Radius parameter: 

10 m < tm < 200 m 

Number of leftover showers: 

322,664 



Figure 3: Measured size distribution of the air showers that have passed the 
quality cuts given in Table [T] The line represents a Gaussian fit to the total 
distribution around the peak, giving a shower size threshold of logj^Q = 
5.92 for the array. 


particles in the Earth’s atmosphere are treated using QGSJET- 
11-04 at high energies and ELUKA ll25l] for energies be¬ 
low 200 GeV. The electromagnetic interactions are treated with 
EGS4 The observation level of the LORA array is set 
to 7.6 m above sea level. Air showers are simulated for pro¬ 
tons and iron nuclei in the energy range of 10^® — 10^® eV, 
assuming a differential energy spectrum with an index —2. The 
showers are weighted to generate a distribution with a spectral 
index —3. Zenith angles are considered in the range 0° — 45°. 
In order to reduce the excessive computing times involved in 
generating the showers, ‘thinning’ is applied at a level of 10“® 
with optimized weight limitation . 


4.2. Detector simulation 


The generated air shower particles are fed into a detector 
simulation code, based on the GEANT4 package ll28ll . which 
allows to calculate the total energy deposition in each detector. 
All properties of the detector, such as the type and the density 
of the scintillator material, the detector geometry as well as the 
effect of the aluminum plates covering the scintillator plates are 
included in the simulation. In order to avoid air showers not 
creating a trigger in the detectors due to the large detector spac¬ 
ing of the LORA array, an additional step is applied to each 


simulated shower before feeding the particles into GEANT4. 
Concentric rings with a radial bin size of 2 m centered around 
the shower axis are constructed, and the total number of par¬ 
ticles contained in each projected ring on the ground is calcu¬ 
lated. All particles in a ring are then distributed uniformly in a 
small square region of area Ag = (1.5 x 1.5) m^ with a LORA 
detector in its center. Depending on the arrival direction of the 
particles, those that hit the detector are allowed to pass through 
GEANT4 and the total energy deposition E’dep in the detector 
is obtained. In the hnal step, the actual amount of energy that 
would have been deposited in the detector is obtained by ap¬ 
plying a correction E’^gp = E^depAg cos 6*/Ar, where 6 is the 
zenith angle of the shower and Ar/ cosO is the projected area 
of the ring on the ground. The somewhat larger area of Ag than 
the actual detector area is used to accommodate particles hitting 
the detector at larger zenith angles. Eor each simulated shower, 
the radial distribution of the energy deposition in the detector, 
averaged over the azimuthal direction in the shower plane, is 
constructed as a function of the distance to the shower axis. 
This method also automatically allows to correct for the effect 
of the shower thinning applied in CORSIKA as the calculation 
takes into account all the particles arriving at the ground. 

Simulations have also been performed to calculate the energy 
deposition of singly charged particles in the detector. Eor that, 
muons of an energy of 4 GeV are considered. Energy depo¬ 
sitions for vertical incident muons and for muons following a 
realistic (observed) arrival direction distribution are obtained. 
The energy deposition distribution for vertical muons gives a 
most probable value of E^vem = 5.3 MeV, while the all-sky 
distribution gives 6.67 MeV. The latter is obtained by also tak¬ 
ing into account a noise level of cr = 1 MeV, which includes a 
contribution from statistical noise, generated by the low number 
of scintillation photons producing a signal and the electronic 
noise. The energy deposition for the all-sky distribution is used 
to calibrate the distribution of the total energy deposition by 
single particles measured with the experiment. Details about 
the calibration are described in Ref. lll9|] . 


4.3. Reconstruction of shower parameters 

Every simulated shower is assigned a random position on the 
ground. The position of the shower axis, and also the detec¬ 
tor coordinates, are then projected in the shower plane. Based 
on the distance of the detector from the position of the shower 
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Figure 4: Combined trigger and reconstruction efficiencies obtained from the simulation of showers induced by protons (left) and iron nuclei (right) as a function of 
the true energy. Different symbols represent different zenith angle bins. See Section l4^ for details. 



Figure 5: Total acceptance for showers induced by protons (squares) and iron 
nuclei (circles) obtained from simulations as function of the true energy. The 
acceptance is calculated for solid angles subtended within 0° — 35°. See Sec- 
tion l4.4l for details. 


axis in the shower plane, the amount of energy deposited in the 
detector is calculated from the radial distribution of energy de¬ 
position given by the simulation. To make the simulation study 
consistent with the analysis of the measured data, the number of 
particles hitting the detectors is obtained in units of VEM (ver¬ 
tical equivalent muons). This is done by first dividing the track- 
length-corrected energy deposition cos 9 by E'vem to ob¬ 
tain the mean number of VEM particles n, hitting the detector. 
To obtain a realistic value, the detector is assigned a number, 
drawn randomly from a Poisson distribution with mean n. This 
last step is necessary to correct for the azimuthal averaging of 
the energy depositions around the shower axis, applied in the 
simulation. The final value rif for the number of VEM particles 
is obtained by adding a random noise, drawn from a Gaussian 
distribution with a standard deviation ct/E^vem- The particle 
density in each detector is obtained by dividing rif by the pro¬ 
jected area of the detector cos 6, where is the actual ge¬ 
ometrical area of the detector. After obtaining the particle den¬ 


sities in the detectors, the reconstruction of air shower parame¬ 
ters is performed similar to the reconstruction of the measured 
air shower data. 

4.4. Trigger and reconstruction efficiencies 

In order to improve the statistics, each simulated shower is 
processed 100 times with the position of the shower axis se¬ 
lected randomly within a circle with a radius of 160 m from 
the center of the array. The fiducial cut of 150 m applied in 
the data analysis is also applied in the calculation of the de¬ 
tector efficiency. A larger radius of 160 m with respect to the 
fiducial cut is necessary to take into account the spillover of 
reconstructed showers across the fiducial boundary due to the 
limited reconstruction accuracy in the position of the shower 
axis which reaches ^ 10 m at a distance of 150 m from the ar¬ 
ray center. Only showers with zenith angles within 0° — 35° are 
considered, and are divided into four different zenith angle bins 
as in the data analysis. Eor each energy and zenith angle bin, 
the trigger efficiency, Ct, is determined by taking the ratio of 
the number of showers that pass through the trigger condition 
listed in Table [T]to the total number of showers generated with 
true shower axis position within the fiducial area. The recon¬ 
struction efficiency, Cr, is calculated as the ratio of the number 
of showers that pass through both the trigger and the quality 
cuts to the total number of triggered showers. Then, the total 
efficiency is obtained as, etot = et Cr- Eigure|4] shows the total 
efficiency for protons (left panel) and iron nuclei (right panel) 
as a function of the true energy for the four zenith angle 
bins; 0° - 15°, 15° - 24°, 24° - 30°, and 30° - 35°. The 
full efficiency of 100% is reached at log]^Q(i?T/GeV) ~ 7.6 
for protons and at ss 7.7 for iron nuclei. 

Eigure |5] shows the total acceptance of the array for primary 
protons and iron nuclei as a function of the true energy. The 
detector acceptance Aacc is defined as the total effective area 
of the array multiplied by the effective viewing angle, and it is 
calculated as, 

A,,,{Et)= [ ^ Ap,oi{9)etot{ET,e,^)dn, (3) 
Jo 
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Figure 6: Left: Comparison of normalized disti'ibutions of radius parameter obtained from the measurements (points) and simulations (thick-solid line: protons 
and thin-dashed line: iron nuclei) for showers with reconstructed sizes log^^Q A^ch > 6.40 and reconstructed zenith angles between 0° and 15°. The inset shows 
Gaussian fits (thin lines) to the distributions around the maximum. Right: Comparison of the averaged lateral distributions between measurements and simulations. 
The measurements (symbols) are the same as shown in Figure [fright panel and the lines (solid: protons and dashed: iron nuclei) are the simulation results for the 
same shower size bin. 


where dfl = sin 0 dO d(j) is the solid angle subtended by an ele¬ 
ment of opening angle between 9 and 9 + d9 and an azimuthal 
width of d(t). flc is the maximum solid angle corresponding to 
the zenith angle cut of 9^ = 35°, Aproj = tt cos 9 is the 
projected geometrical area of the array at an inclination 9 with 
Rc = 150 m representing the fiducial radial cut applied in the 
analysis, and the total efficiency, etot. is given as a function of 
{Et,9, (/)). Assuming azimuthal symmetry of et, the integral in 
Equation[3is discretized in zenith angle bins and can be rewrit¬ 
ten as, 

AacciEr) = —(cos20k — cos20k+i), 

(4) 

where k denotes the zenith angle bins, ng = 4 is the number 
of zenith angle bins considered, and 0k and 0k-i-i represent the 
low-bin and high-bin edges of each zenith angle bin respec¬ 
tively. 

5. Comparison between simulations and measurements 

In Figure |6] (left panel), the normalized distribution of radius 
parameters for the simulated showers with reconstructed sizes 
logiQ A"ch > 6.40 is compared with the measurements for the 
zenith angle range of 0° — 15°. The points in the figure repre¬ 
sent the measurements and they are the same as shown in Fig- 
ure|2](left panel). The distribution for iron nuclei (thick-dashed 
line) shows a systematic shift towards larger ru with respect to 
the proton showers (thick-solid line) which is expected due to 
the difference in the shower development between proton and 
iron primaries. Showers induced by iron nuclei are generated 
higher up in the atmosphere, resulting in a larger spread (which 
implies larger tm values) on the ground, relative to the pro¬ 
ton showers. Although both the simulated distributions follow 
a similar shape as the measured distribution, they are not in 


full agreement with the data. But, overall, the proton distribu¬ 
tion seems to be relatively closer to the data. The inset shows 
a closer view for the region around the maximum between 12 
and 48 m. The lines represent fits to the distributions using a 
Gaussian function. The Gaussian peaks for the simulated dis¬ 
tributions obtained from the fits are (30.51 ± 0.01) m for the 
proton distribution and (34.38 ± 0.02) m for the iron distribu¬ 
tion. The value for the proton distribution is found to be quite 
close to the peak value of (30.33 ± 0.13) m obtained for the 
data. 

Figure |6] (right panel) shows a comparison of the averaged 
lateral distribution between simulations and measurements for 
a reconstructed shower size in the range of 6.40 < logj^Q A^ch < 
7.75. The measurements (points) are the same as already shown 
in Figure |2] (right panel). The iron distributions (dashed lines) 
are found to be slightly flatter than the proton distributions 
(solid lines), which is expected due to the larger tm values for 
iron showers as explained above. Although both the simulated 
proton and iron distributions are consistent with the data within 
the experimental uncertainties, the proton distributions seem to 
agree better as the iron distributions tend to show some sys¬ 
tematic deviation from the measurements above a distance of 
70 m from the shower axis. A test of the comparison 
between the simulations and the measurements gives reduced 
values within the range of ~ 1.15 — 1.26 for protons and 
~ 1.33 — 1.86 for the case of iron nuclei. The better agreement 
of the measurements with the proton distributions is expected 
because the air shower particles measured by FORA are mostly 
dominated by electrons rather than muons, which makes the 
measurements biased towards protons. 

6. Energy calibration 

The measured shower size can be converted into the energy 
of the primary particle using a conversion relation obtained 
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Figure 7: Two-dimensional histogram for the reconstructed shower size A^ch true energy E^, obtained from simulations for showers induced by protons (left) 
and iron nuclei (right) for a zenith angle bin of 0° — 15°. The distributions are weighted to an energy spectrum of index -3. Each point represents the peak in 
the true energy distribution for a logj^g A^ch bin width of 0.15 (see Section[^for details). The lines represent fits to the points using Equation[6|in the range of 
6.4 < logj^gA^ch < 8-3 for protons and 6.3 < log^gA^ch < 8.3 for iron nuclei. The fit parameters are listed in Table 2. 
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Figure 8: True energy distribution for the reconstructed size bin of 6.70 < logjQ77t.h < 6.85 for protons (left panel) and iron nuclei (right panel) for the zenith 
range 0° — 15°. The distributions are weighted to an energy spectrum of index -3. The lines represent fits using a skewed Gaussian function given by Equation|^ 
which involves four parameters (see Section|^for details). Main fit parameters (pi, P 21 P 3 ) are shown. 


from simulations. Simulated showers are stored in a two dimen¬ 
sional log-log histogram in reconstructed size and true energy. 
Such a histogram is shown in Figure [T] for showers induced by 
protons (left panel) and iron nuclei (right panel) for the zenith 
angle range of 0° — 15°. The color profile represents the weight 
of the distribution. The distribution is broader for proton show¬ 
ers, which is mainly due to the large intrinsic fluctuations of 
proton showers. From simulations, the fluctuations in the true 
shower size for proton showers within the 0° — 15° zenith an¬ 
gle bin are found to be ^ (30 — 45)% while the uncertainty 
due to the reconstruction is in the range of ^ (12 — 19)% for 
the energy region of our interest. For iron induced showers, the 
intrinsic size fluctuation is only ^ (17 — 22)%, while the re¬ 
construction accuracy remains almost the same as that of the 
proton induced showers. Another major difference is that for 
the same reconstructed shower size, iron showers have higher 
energies than the protons. This is related to the shallower pene¬ 
tration depth of iron induced showers in the atmosphere, which 


leads to an increased attenuation of electrons before they can 
reach the ground. 

The distributions in Figure|2]are binned in A^ch. taking a log¬ 
arithmic bin size of 0.15, and profile plots of the true energy 
as function of Nch are generated. The profile plots are repre¬ 
sented by the solid points in Figure [T] Each point in the plots 
represents the peak of the energy distribution for each size bin, 
and the uncertainty on each point corresponds to the spread 
of the energy distribution which is described in detail in the 
following. Figure 0 shows the energy distribution for the bin 
slice of logj^o -^ch = 6.70 — 6.85 for both, showers induced 
by protons (left panel) and iron nuclei (right panel). The dis¬ 
tribution for protons is not symmetric about its mean and is 
found to be more extended to lower energies. This can be un¬ 
derstood as more contamination from low-energy showers in a 
given size bin than from higher energies which is caused by the 
larger intrinsic fluctuations of low-energy showers. The level 
of contamination depends on the assumed slope of the primary 
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Table 2: Fit parameters for showers induced by protons and iron nuclei obtained by fitting Equation|6]to the size-energy profile plots for different zenith angle bins 
between 0° and 35°. A slope 7 s = —3 of the differential cosmic-ray energy spectrum has been adopted in the simulations. 


Zenith angle 



Protons 





Iron nuclei 
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a 



b 



a 



b 


0° ■ 

- 15° 

0.980 

± 

0.683 

0.922 

± 

0.089 

1.747 

± 

0.361 

0.853 

± 

0.048 

15° 

-24° 

1.234 

± 

0.766 

0.898 

± 

0.099 

1.801 

± 

0.370 

0.858 

± 

0.049 

24° 

-30° 

1.315 

± 

0.815 

0.901 

± 

0.101 

1.726 

± 

0.319 

0.885 

± 

0.041 

O 

O 

CO 

-35° 

1.667 

± 

0.692 

0.873 

± 

0.091 

1.982 

± 

0.366 

0.866 

± 

0.048 


cosmic-ray spectrum in the simulation. The peaks of the distri¬ 
butions are obtained by fitting with a skewed Gaussian function. 
The skewed Gaussian distribution function used in the present 
analysis is given by, 


f{x) = — exp 
Pi 


( {X-P2f \ 

\ ) 


1 + erf 


[ P3ix-P2)\ 

(5) 


where x = logj^Q (i^x/GeV), pq is the normalisation constant, 
Pi and p2 represent measures of the spread and the position of 
the distribution respectively, and p^ is the skewness parameter 
of the distribution function. 

The lines in Figure [^represent the fitted functions. The im¬ 
portant fit parameters (pi, P 2 , P 3 ) are also shown. For the distri¬ 
bution of iron-induced showers, it can be noticed that the value 
of p 3 is close to zero, indicating that the distribution closely re¬ 
sembles a normal Gaussian distribution. The uncertainties in 
the profile plots shown in Figure |7] are obtained by taking the 
difference between the energies corresponding to the full width 
at half maximum (FWHM) and the peak energy of the energy 
distribution for each size bin. The uncertainties obtained are 
asymmetric for the proton distribution while for iron induced 
showers, they are almost symmetric. The iVch values of the 
profile plots shown in Figure [T] are calculated as the weighted 
mean of the size distribution within each size bin. These size 
values are found to be slightly smaller than the bin centers. 

To obtain the size-energy relation, each profile plot is fitted 
using the following function. 


logio Et = a + b logio ^ch (6) 

where Et denotes the true energy, and a and b are the fit pa¬ 
rameters. The fit is performed only in the size region where a 
reliable fit of the true energy distribution, as shown in Figure 
[8] could be performed. This corxesponds to a size region of 
logj^o -^ch = 6.4 — 8.3 for both the type of particles. The pro¬ 
file plots as well as the fitted functions for all the zenith angle 
ranges are shown in Figure 0 for protons (left panel) and iron 
nuclei (right panel). 

From these figures, it can be noticed that for the same shower 
size, primary energies at larger zenith angles are larger than at 
smaller angles. In other words, it requires a higher energy at 
larger zenith angles to generate the same number of particles 
on the ground as at lower zenith angles. This is due to higher 
attenuation of air shower particles at larger zenith angles as the 
showers pass through a longer column depth of air in the atmo¬ 
sphere. The values of the a and b parameters obtained from the 


fits for the four zenith angle ranges are listed in Table |2l Using 
these values, for any simulated or measured shower for which 
the reconstructed arrival direction and the reconstructed size are 
known, the primary cosmic-ray energy can be reconstructed us¬ 
ing the relation 


logio EYi = a + b logio ^ch, (7) 

where E-^ denotes the reconstructed energy. 

7. Energy resolution and systematic uncertainties 

In this section, we present details about the accuracy of the 
reconstructed energies and the uncertainties that have to be 
considered for the reconstruction of the cosmic-ray intensity. 
The accuracy depends on the variation of the true shower size 
caused by the intrinsic shower-to-shower fluctuations in the at¬ 
mosphere and also on the accuracy in the reconstruction of the 
shower size. 

7.7. Energy resolution 

For each size bin in the Et — Nch profile plot, reconstructed 
energies (7 ?r) are obtained for every simulated shower, and a 
distribution of the differences between the true energies and the 
reconstructed energies (T^x — 77 r) is generated. The distribu¬ 
tion obtained is similar to the one shown in Figure [S] except 
for a shift in the peak position to the left by an interval equal 
to the value of the reconstructed energy. The peak position and 
the spread of these distributions are obtained correspondingly. 
The peak represents the systematic uncertainty due to energy 
calibration, while the spread corresponds to the energy reso¬ 
lution. Their values expressed as fraction of the reconstructed 
energies are shown in Figure [10] as function of the shower size 
for the zenith angle range of 0 = 0° — 15°. The resolution 
is in the range of ^ 28% — 48% for proton induced showers 
and ~ 12% — 32% for showers induced by iron nuclei. The 
systematics are within 12% for protons and within 8% for iron 
nuclei. At 6* = 30° — 35°, the uncertainty in energy for protons 
increases to the range of ^ 37% — 65% in resolution and to 
^ 20% in systematics. For iron nuclei, the uncertainty remains 
almost the same up to 6* = 35°. 

7.2. Systematic uncertainty in energy 

The systematic uncertainty in energy shown in Figure [T0| is 
associated with the energy calibration performed using Equa¬ 
tion [T] Other main sources of systematic uncertainty in energy 
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Figure 9: Reconstructed size and true energy (-Ex) relation for proton (left) and iron (right) showers for all zenith angle ranges up to 35°. 
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Figure 10: Accuracy in the reconstructed energies (Ej^) for showers induced 
by protons (squares) and iron nuclei (circles) as a function of the reconstructed 
size Ech for the zenith angle bin of 0° — 15°. The filled points represent 
the energy resolution and the the empty points are the systematic uncertainties 
resulting from the energy calibration. See Section|TT]for details. 


include the assumed slope of the primary cosmic-ray spectrum 
in the CORSIKA simulation, the VEM peak obtained from the 
detector simulation and the hadronic interaction models. Thus, 
the calibration parameters, listed in Table |2] also depend on the 
choice of simulation parameters. 

A part of the systematic uncertainties are obtained by chang¬ 
ing the values of the slope and the VEM peak in the simulations 
within reasonable limits, and by comparing the newly recon¬ 
structed energies with the energies obtained using the fixed pa¬ 
rameters given in Table|2] Eor the slope of the energy spectrum, 
simulated showers with an original slope 7 s = —2 are weighted 
to generate distributions for 7 s = —2.5 and 7 s = —3.5. Then, 
following the same procedure as described in Section en¬ 
ergy calibrations are performed separately for the two different 
slopes and calibration parameters are obtained. The differences 
between the energies reconstructed with the new parameters 
and the ones reconstructed using the parameters given in Ta¬ 
ble |2]gives the systematic uncertainty due to the spectral slope. 


The uncertainties are found to be within (-1-6%, —9%) for pro¬ 
tons and within ±2% for iron nuclei. 

Erom the detector simulation, it has been observed that 
adding noise to the deposited energy in the detector at the level 
of 1 MeV (see Section |4^ leads to around 10% positive shift 
in the value of the most probable energy deposition Evem in 
the detector for vertical incident muons. The 10% increase in 
Evem will lead to a decrease in the shower size and subse¬ 
quently to an increase in the reconstructed energy by ~ 10%. 
The average systematic shift in the reconstructed energy due to 
this uncertainty in VEM calibration is obtained to be ~ -f 10% 
for showers induced by either protons or iron nuclei. The dif¬ 
ferent systematic uncertainties obtained are shown in EigurefTTI 
as a function of the reconstructed energy for showers induced 
by protons (left panel) and iron nuclei (right panel). Eor proton 
showers, the total systematic uncertainty, obtained by adding 
the individual systematic components in quadrature, is found 
to be within ~ (-1-20%, —10%) and for iron showers, the to¬ 
tal systematic is within (-1-10%,—5%). At larger zenith an¬ 
gles, the total systematic for protons increases slightly, reaching 
~ (-1-22%, —15%) at 0 = 30° — 35°, while for iron nuclei, the 
uncertainty remains almost unchanged. 

7.3. Systematic uncertainty in intensity 

Any systematic uncertainty in energy results in a systematic 
shift in the reconstructed cosmic-ray flux intensity. To estimate 
the systematic uncertainty in intensity due to the energy calibra¬ 
tion, the reconstructed energies are determined using Equation 
|2]for all simulated showers with 7 s = —3 that pass through all 
selection and quality cuts as listed in Table [T] The distribution 
of the reconstructed energies is compared to the distribution of 
the true energies, and the systematic uncertainty in intensity is 
calculated as (It — Ik)/Ik for each energy bin, where /x and 
Ik represent the number of showers per bin in the true and re¬ 
constructed energy distributions respectively. 

Eor the systematic effect due to the uncertainties in the spec¬ 
tral slope and the VEM calibration, the energy calibration de¬ 
termined in their respective cases are applied to the simulated 
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Figure 11: Systematic uncertainties in the reconstructed energy (Er) for showers induced by protons (left panel) and iron nuclei (right) for the zenith angle bin 
of 0° — 15° as a function of i?R. The systematic uncertainties due to the energy calibration (thick solid lines) are the same as shown in Figure [Tol but plotted as 
function of Er. They are calculated using the parameters set given in Table[2]for 7 a = —3. The blue band represents the uncertainty resulting from changing the 
specttal slope from —2.5 to —3.5. The dashed line is due to the uncertainty involved in the VEM calibration and the shaded-striped region is the total uncertainty. 
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Figure 12: Systematic uncertainties in intensity obtained from simulations of showers induced by protons (left panel) and iron nuclei (right panel) for the zenith 
angle bin of 0° — 15°. All lines/bands (except the dotted band) have the same representation as in Figure [171 The dotted band represents the uncertainty expected 
due to the hadronic interaction model which is taken as 12.5%. 


showers for 7 s = —3 and the distributions of the newly re¬ 
constructed energies are compared with the old distribution ob¬ 
tained using the parameters given in Table |2] 

Figure [12] shows the different systematic uncertainties in in¬ 
tensity that have been obtained for protons (left panel) and iron 
nuclei (right panel). The thick solid line represents the sys¬ 
tematic uncertainty due to the energy calibration, the blue band 
represents the contribution due to the spectral slope, the dashed 
line is the VEM contribution, and the shaded-striped region 
represents the total systematic uncertainty. For energies above 
log]^Q(£’R/GeV) ~ 7.2, the systematic uncertainty due to the 
energy calibration is found to be within ~ (-1-30%,—10%) 
for proton showers and within ^ (-1-20%,—10%) for show¬ 
ers induced by iron nuclei. The systematic uncertainty due 
to the spectral slope is within ^ (-1-40%,—15%) for pro¬ 
tons, and within ^ (-1-12%, —18%) for iron nuclei except at 
log]^Q(i7/GeV) ^ 8.6 where the uncertainty reaches ~ 30%. 
The systematic uncertainty associated with the VEM calibra¬ 


tion is found to be within -|-30% for both types of nuclei. A con¬ 
tribution of 12.5% due to the uncertainty in the hadronic inter¬ 
action model 1 3^ 31] is also included in Figure[T2] For protons, 
the total systematic uncertainty above log]^g(i7R/GeV) ~ 7.2 
is within ~ (-1-60%, —25%), and for iron nuclei, the total un¬ 
certainty is within ~ (-1-38%, —20%). 


8. Measured cosmic-ray energy spectrum 

For all high-quality LORA data, reconstructed energies are 
determined on shower-by-shower basis, and a distribution of 
reconstructed energies is built taking a logarithmic bin size of 
0.15. From the distribution, the differential cosmic-ray spec¬ 
trum {dl/dE) is obtained by folding in the total acceptance of 
the LORA array Aacc (Figure |5|l and the total observation time 
Tobs as follows. 

An 

Xe 


1 


^accTobs 


(8) 
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Figure 13: Left: All-particle cosmic-ray energy spectrum measured with LORA, assuming that cosmic rays are only protons (squares) and iron nuclei (filled 
circles). The error bars represent statistical uncertainties and the shaded areas represent systematic uncertainties. The lines represent single power law fits to the 
measurements, excluding the highest three energy bins. Right: LORA measurements compared to the all-particle energy spectrum from IceTop (crosses) and 
KASCADE-Grande (empty circles) measurements. 


where the subscript i denotes the i*** energy bin and An is the 
number of showers in an energy bin of width Ai?. For con¬ 
structing the spectrum, only the energy region that has a total 
(trigger and reconstruction) efficiency greater than 98% is used. 
This corresponds to an energy of 1.9 x 10^ GeV for protons and 
2.7 X 10^ GeV for iron nuclei (see Figure|4]l. 

Figure [TS] (left panel) shows the reconstructed energy spec¬ 
trum multiplied by E^, assuming that cosmic rays are only pro¬ 
tons or iron nuclei. The spectrum is given in the energy range 
of (1.9 X 10^ — 1.2 X 10®) GeV for protons, and in the range 
of (2.7 X 10^ — 1.7 X 10®) GeV for iron nuclei. The mea¬ 
sured values along with the uncertainties are listed in Table [3] 
The measured spectra cannot be described by single power laws 
over the full energy range because of the structures present in 
the spectra, particularly the dip at ^ 6 x 10® GeV. A power 
law fit to the measured spectra data below 5 x 10® GeV gives 
spectral index values of 7 p = —3.18 ± 0.13 for protons and 
7Fe = —3.22 ± 0.08 for iron nuclei. 

In Figure [13] (right panel), our measured spectra are com¬ 
pared with the all-particle spectra measured with the IceTop 
i29l] and KASCADE-Grande experiments. Both their 
spectra lie between our reconstructed spectra, which is expected 
in the case of a mixed cosmic-ray composition. They are close 
to our proton spectrum at ~ 2 x 10^ GeV, and become closer to 
our iron spectrum as the energy increases. This might be an in¬ 
dication of a change in the mass composition of cosmic rays in 
the energy region between 10^ and 10® GeV, which is expected 
as due to a transition from a Galactic to an extragalactic origin 
of cosmic rays. 

9. Conclusion and outlook 

We have conducted a detailed energy reconstruction study for 
the extensive air showers measured with the LORA particle de¬ 
tector array. Important parameters such as the energy resolution 
of the array and the systematic uncertainty of the reconstructed 
energy have been obtained. The energy resolution is found to be 


in the range of ^ 28 — 48% for showers induced by protons and 
^ 12 — 32% for iron nuclei. The total systematic uncertainty of 
the reconstructed energy is within ~ (-1-20%, —10%) for pro¬ 
tons and within ~ (-1-10%, —5%) for iron nuclei. Applying 
the reconstruction method to the measured data, the all-particle 
cosmic-ray energy spectrum has been obtained, assuming that 
cosmic rays are only constituted of protons or iron nuclei for 
energies above ~ 10^® eV with a systematic uncertainty in in¬ 
tensity of ^ 20 — 60%. Our future effort will concentrate on 
combining the energy measurement of LORA with the compo¬ 
sition measurement from the LOLAR radio antennas to deter¬ 
mine an all-particle energy spectrum, taking into account the 
actual cosmic-ray composition. 

Especially the primary energy determined using the energy 
calibration given here is being used in the reconstruction of air 
shower properties with the radio data from LOLAR. Calcula¬ 
tion of energy calibration parameters for higher zenith angles 
above ~ 40° is underway. This is particularly important for 
the LOLAR radio measurements where a significant fraction of 
showers have been observed at larger zenith angles. At present, 
the small size of the LORA array effectively limits the effective 
area of LOLAR. Efforts are ongoing to expand the size of the 
array to exploit the full potential of LOLAR. 
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Table 3: Values of the measured cosmic-ray spectrum, assuming that cosmic rays are only protons or iron nuclei. The energies are given in GeV and the intensities 
along with the statistical and the systematic uncertainties are given in units of l/(m^ sr s GeV). 


Energy Intensity ± stat. ± sys. uncertainties [l/(m^ sr s GeV)] 
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